#Importando bibliotecas

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.0     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.1     ✔ tibble    3.2.0
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(readxl)
library(dplyr)
library(clipr)
## Welcome to clipr. See ?write_clip for advisories on writing to the clipboard in R.
library(did)
library(ggplot2)
library(modelsummary)
## `modelsummary` 2.0.0 now uses `tinytable` as its default table-drawing
##   backend. Learn more at: https://vincentarelbundock.github.io/tinytable/
## 
## Revert to `kableExtra` for one session:
## 
##   options(modelsummary_factory_default = 'kableExtra')
## 
## Change the default backend persistently:
## 
##   config_modelsummary(factory_default = 'gt')
## 
## Silence this message forever:
## 
##   config_modelsummary(startup_message = FALSE)
setwd("C:/Users/Ana Flávia/OneDrive - Insper - Institudo de Ensino e Pesquisa/Economia/Quinto Semestre/Microeconomia IV/Trabalho Final")

1 Ajustando as bases de dados

1.1 Base: mortes por overdose por estado

mortes = read.csv("NCHS_-_Drug_Poisoning_Mortality_by_State__United_States_20240417 (1).csv")

mortes = mortes %>%
  select(State, Year, Sex, Age.Group, Race.and.Hispanic.Origin, Deaths, Population) %>% 
  mutate(Year = as.character(Year)) %>% 
  filter(Sex == "Both Sexes", Age.Group == "All Ages", 
         Race.and.Hispanic.Origin == "All Races-All Origins") %>% 
  select(State, Year, Deaths, Population) %>% 
  mutate(Overdoses = Deaths) %>% 
  select(State, Year, Overdoses, Population)%>% 
  rename(year = Year) 

mortes <- mortes %>%
  filter(State != "United States",
         year != "1999",
         year != "2000",
         year != "2001",
         year != "2002")

view(mortes)

1.2 Base: IDH por estado

IDH = read.csv("GDL-Subnational-HDI-data.csv", sep = ";", quote = "", 
               header = TRUE)

IDH = IDH %>% 
  rename(State = Region) %>% 
  select(State, X2003, X2004, X2005, X2006,X2007, X2008, X2009, 
         X2010, X2011, X2012, X2013, X2014, X2015, X2016) %>% 
  slice(-1) %>% 
  pivot_longer(!State, names_to = "year", values_to = "IDH") %>% 
  mutate(year = gsub("X","", year))

view(IDH)

1.3 Base: desemprego por estado

unemp = read_xls("emp-unemployment.xls", sheet = "States")
## New names:
## • `` -> `...2`
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`
## • `` -> `...15`
## • `` -> `...16`
## • `` -> `...17`
## • `` -> `...18`
## • `` -> `...19`
## • `` -> `...20`
## • `` -> `...21`
## • `` -> `...22`
## • `` -> `...23`
## • `` -> `...24`
## • `` -> `...25`
## • `` -> `...26`
## • `` -> `...27`
## • `` -> `...28`
## • `` -> `...29`
## • `` -> `...30`
## • `` -> `...31`
## • `` -> `...32`
## • `` -> `...33`
## • `` -> `...34`
## • `` -> `...35`
## • `` -> `...36`
## • `` -> `...37`
## • `` -> `...38`
## • `` -> `...39`
## • `` -> `...40`
## • `` -> `...41`
unemp = unemp %>% 
  slice(-(1:4)) %>% 
  slice(-(54:62))

header = unemp %>% 
  slice(1) %>% 
  unlist()

unemp = unemp %>% 
  slice(-1) %>% 
  rename_with(~header) %>% 
  slice(-1) %>% 
  rename_with(
    ~ paste0("X", .)) %>% 
  rename(State = XArea,
         Fips = XFips) %>% 
  select(State, X2003, X2004, X2005, X2006, X2007, X2008, 
         X2009, X2010, X2011, X2012, X2013, X2014, X2015, X2016) %>% 
  pivot_longer(!State, names_to = "year", values_to = "Desemprego") %>% 
  mutate(year = gsub("X","", year))
    
view(unemp)

1.4 Base: renda por estado

income = read_xlsx("h08.xlsx")
## New names:
## • `` -> `...2`
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`
## • `` -> `...15`
## • `` -> `...16`
## • `` -> `...17`
## • `` -> `...18`
## • `` -> `...19`
## • `` -> `...20`
## • `` -> `...21`
## • `` -> `...22`
## • `` -> `...23`
## • `` -> `...24`
## • `` -> `...25`
## • `` -> `...26`
## • `` -> `...27`
## • `` -> `...28`
## • `` -> `...29`
## • `` -> `...30`
## • `` -> `...31`
## • `` -> `...32`
## • `` -> `...33`
## • `` -> `...34`
## • `` -> `...35`
## • `` -> `...36`
## • `` -> `...37`
## • `` -> `...38`
## • `` -> `...39`
## • `` -> `...40`
income = income %>% 
  slice(-(1:60)) 

header = income %>% 
  slice(1) %>% 
  unlist()

income = income %>% 
  slice(-1) %>% 
  rename_with(~header) %>% 
  slice(-1) %>% 
  rename_with(
    ~ paste0("X", .)) %>% 
  rename(State = XState) %>% 
  select(State,X2003, X2004,X2005,X2006, X2007, X2008, X2009, 
         X2010, X2011, X2012, X2013, X2014, X2015, X2016) %>% 
  pivot_longer(!State, names_to = "year", values_to = "Income") %>% 
  mutate(year = gsub("X","", year))

view(income)

1.5 Base: leis maconha por estado

leis  = read_xlsx("BASE_LEISMACONHA.xlsx", sheet = "RML")

leis = leis %>% 
  select(-(2:32)) %>% 
  select(-(16:23))

leis = leis %>% 
  rename(State = "Estados/Anos") %>% 
  pivot_longer(!State, names_to = "year", values_to = "Tratamento") %>% 
  mutate(Tratamento = as.numeric(Tratamento))

view(leis)

1.6 Base: uso de drogas (maconha e cocaína) por estado

drug_users = read_xlsx("Drug users.xlsx")

drug_users <- drug_users %>% 
  mutate(year = as.character(year))


view(drug_users)

1.7 Base unificada

base_TF = left_join(mortes, IDH, by = c("State", "year"))
base_TF = left_join(base_TF, unemp, by = c("State", "year"))
base_TF = left_join(base_TF, income, by = c("State", "year"))
base_TF = left_join(base_TF, leis, by = c("State", "year"))
base_TF = left_join(base_TF, drug_users, by = c("State", "year"))

base_TF <- base_TF %>% 
  mutate(year = as.numeric(year),
         weed_users = as.numeric(weed_users))

base_TF <- base_TF %>% 
  mutate(Fips = case_when(State == 'Alabama' ~ '01000',
                          State == 'Alaska' ~'02000',
                          State == 'Arizona' ~'04000',
                          State == 'Arkansas' ~'05000',
                          State == 'California' ~'06000',
                          State == 'Colorado' ~'08000',
                          State == 'Connecticut' ~'09000',
                          State == 'Delaware' ~'10000',
                          State == 'District of Columbia' ~'11000',
                          State == 'Florida' ~'12000',
                          State == 'Georgia' ~'13000',
                          State == 'Hawaii' ~'15000',
                          State == 'Idaho' ~'16000',
                          State == 'Illinois' ~'17000',
                          State == 'Indiana' ~'18000',
                          State == 'Iowa' ~'19000',
                          State == 'Kansas' ~'20000',
                          State == 'Kentucky' ~'21000',
                          State == 'Louisiana' ~'22000',
                          State == 'Maine' ~'23000',
                          State == 'Maryland' ~'24000',
                          State == 'Massachusetts' ~'25000',
                          State == 'Michigan' ~'26000',
                          State == 'Minnesota' ~'27000',
                          State == 'Mississippi' ~'28000',
                          State == 'Missouri' ~'29000',
                          State == 'Montana' ~'30000',
                          State == 'Nebraska' ~'31000',
                          State == 'Nevada' ~'32000',
                          State == 'New Hampshire' ~'33000',
                          State == 'New Jersey' ~'34000',
                          State == 'New Mexico' ~'35000',
                          State == 'New York' ~'36000',
                          State == 'North Carolina' ~'37000',
                          State == 'North Dakota' ~'38000',
                          State == 'Ohio' ~'39000',
                          State == 'Oklahoma' ~'40000',
                          State == 'Oregon' ~'41000',
                          State == 'Pennsylvania' ~'42000',
                          State == 'Rhode Island' ~'44000',
                          State == 'South Carolina' ~'45000',
                          State == 'South Dakota' ~'46000',
                          State == 'Tennessee' ~'47000',
                          State == 'Texas' ~'48000',
                          State == 'Utah' ~'49000',
                          State == 'Vermont' ~'50000',
                          State == 'Virginia' ~'51000',
                          State == 'Washington' ~'53000',
                          State == 'West Virginia' ~'54000',
                          State == 'Wisconsin' ~'55000',
                          State == 'Wyoming' ~'56000')) %>% 
           mutate(Fips = as.numeric(Fips))

base_TF<-base_TF %>% 
  mutate(taxa_overdose = (Overdoses/Population)*1000,
         taxa_weed = (weed_users/Population)*1000,
         taxa_cocain = (cocain_users/Population)*1000)

view(base_TF)

2 Descritivas

2.1 Tabelas

#pré tratamento (2012)
base_TF %>% 
  filter(State %in% c("Colorado", "Nevada")) %>%
  filter(year %in% c("2003":"2011")) %>%
  rename(Cocaina = cocain_users,
         Maconha = weed_users,
         Renda = Income, 
         Habitantes = Population) %>% 
  select(State, Overdoses, Habitantes, IDH, Desemprego, Renda,
         Cocaina, Maconha) %>% 
  datasummary(formula = ~ Overdoses + Habitantes + IDH 
                   + Desemprego + Renda + Cocaina + Maconha
                ~ State*(Mean + SD))
tinytable_rz64fbm498i8wovpvtqc
Colorado Nevada
Mean SD Mean SD
Overdoses 667.11 115.80 491.67 93.79
Habitantes 4807537.89 208834.96 2545901.44 169886.65
IDH 0.93 0.00 0.90 0.00
Desemprego 5.97 1.79 7.39 4.03
Renda 76971.11 3454.54 69088.89 4295.52
Cocaina 152024.49 19155.18 55942.56 7872.71
Maconha 713874.37 122214.88 302510.13 38672.94
#pré tratamento (2014)
base_TF %>% 
  filter(State %in% c("Alaska", "Oregon")) %>%
  filter(year %in% c("2003":"2013")) %>% 
  rename(Cocaina = cocain_users,
         Maconha = weed_users,
         Renda = Income, 
         Habitantes = Population) %>% 
  select(State, Overdoses, Habitantes, IDH, Desemprego, Renda,
         Cocaina, Maconha) %>% 
  datasummary(formula = ~ Overdoses + Habitantes + IDH 
                   + Desemprego + Renda + Cocaina + Maconha
                ~ State*(Mean + SD))
tinytable_uvrhds5jpr7vkienbhzv
Alaska Oregon
Mean SD Mean SD
Overdoses 98.00 22.66 456.91 61.25
Habitantes 692375.27 29618.52 3748458.18 133901.36
IDH 0.93 0.01 0.91 0.01
Desemprego 7.19 0.54 7.88 2.04
Renda 79466.36 3991.66 65460.91 2490.47
Cocaina 16466.17 3102.24 83763.70 15332.42
Maconha 121482.87 18702.32 594795.02 100268.74
#summary geral
base_TF %>% 
  filter(!State %in% c("Alaska", "Oregon","Colorado", "Nevada")) %>% 
  rename(Cocaina = cocain_users,
         Maconha = weed_users,
         Renda = Income, 
         Habitantes = Population) %>% 
  select(State, Overdoses, Habitantes, IDH, Desemprego, Renda,
         Cocaina, Maconha) %>% 
  datasummary(formula = ~ Overdoses + Habitantes + IDH 
                   + Desemprego + Renda + Cocaina + Maconha
                ~ Mean + SD)
tinytable_5c0ckozb848nu5ojx9x8
Mean SD
Overdoses 804.42 855.61
Habitantes 6280125.08 6965799.14
IDH 0.91 0.02
Desemprego 5.99 2.00
Renda 66166.26 10341.97
Cocaina 125089.51 151932.98
Maconha 699701.09 842269.89

2.2 Gráficos

#gráficos com os tratamentos
base_TF %>% 
  filter(State %in% c("Alaska", "Colorado", "Oregon",
                      "Washington")) %>% 
  filter(year %in% c("2003":"2013")) %>% 
  ggplot(aes(x = year, y = taxa_overdose))+
  geom_line(color = "#336600")+
  facet_wrap(~State) +
  labs(title = "Overdose para os tratados", 
       x = "Ano", y = "Taxa de overdose por 1000 habitantes")

base_TF %>% 
  filter(State %in% c("Alaska", "Colorado", "Oregon",
                      "Washington")) %>% 
  filter(year %in% c("2003":"2013")) %>% 
  ggplot(aes(x = year, y = taxa_weed))+
  geom_line(color = "#336600")+
  facet_wrap(~State) +
  labs(title = "Usuários de maconha nos grupos tratados", 
       x = "Ano", y = "Usuários de maconha por 1000 habitantes")

base_TF %>% 
  filter(State %in% c("Alaska", "Colorado", "Oregon",
                      "Washington")) %>% 
  filter(year %in% c("2003":"2013")) %>% 
  ggplot(aes(x = year, y = taxa_cocain))+
  geom_line(color = "#336600")+
  facet_wrap(~State) +
  labs(title = "Usuários de cocaína nos grupos tratados", 
       x = "Ano", y = "Usuários de cocaína por 1000 habitantes")

#gráfico de geral
medias <- base_TF %>% 
  filter(!State %in% c("Alaska", "Oregon","Colorado", "Nevada")) %>% 
  group_by(year) %>% 
  summarise(medias_overdose = mean(taxa_overdose),
            medias_maconha = mean(taxa_weed),
            medias_cocaina = mean(taxa_cocain))

medias %>% 
  ggplot(aes(x = year, y = medias_overdose))+
  geom_line(color = "#336600") +
  labs(title = "Médias de overdose nos EUA - geral",   
       x = "Ano", y = "Médias de Overdose por 1000 habitantes")

medias %>% 
  ggplot(aes(x = year, y = medias_maconha))+
  geom_line(color = "#336600") +
  labs(title = "Médias de usuários de maconha nos EUA - geral",   
       x = "Ano", y = "Médias de usuários de maconha por 1000 habitantes")

medias %>% 
  ggplot(aes(x = year, y = medias_cocaina))+
  geom_line(color = "#336600") +
  labs(title = "Médias de usuários de cocaína nos EUA - geral",   
       x = "Ano", y = "Médias dos usuários de cocaína por 1000 habitantes")

3 Regressão

3.1 Regressão de morte por overdoeses

#tentando rodar o sttagerd
primeiros_tratamentos <- base_TF %>%
  filter(Tratamento > 0) %>%
  group_by(State) %>%
  summarise(primeiros_tratamentos = min(year))

view(primeiros_tratamentos)

base_TF <- base_TF %>%
  mutate(year = as.numeric(year),
         Tratamento = as.numeric(Tratamento),
         Fips = as.numeric(Fips)) %>% 
  left_join(primeiros_tratamentos, by = "State") %>%
  mutate(Tratado = ifelse(Tratamento > 0, primeiros_tratamentos, 0)) %>% 
  arrange(State, year) %>%
  fill(Tratado, .direction = "down") %>% 
  select(-primeiros_tratamentos)

base_TF <- base_TF %>% 
  mutate(Tratado = case_when(State == "Alaska"~"2014",
                             State == "California"~"2016",
                             State == "Colorado"~"2012",
                             State == "Maine"~"2016",
                             State == "Massachusetts"~"2016",
                             State == "Nevada"~"2016",
                             State == "Oregon"~"2014",
                             State == "Washington"~"2012",
                             .default = "0")) %>% 
  mutate(Tratado = as.numeric(Tratado))

view(base_TF)

#filtrei esses 4 estados pq não tem pós tratamento pra eles
base_2 <- base_TF %>% 
  filter(State != "California",
         State != "Maine",
         State != "Massachusetts",
         State != "Nevada")
view(base_2)

overdose_death <- att_gt(yname = "taxa_overdose",
                   tname = "year",
                   idname = "Fips",
                   gname = "Tratado",
                   data = base_2)
## Warning in pre_process_did(yname = yname, tname = tname, idname = idname, : Be aware that there are some small groups in your dataset.
##   Check groups: 2012,2014.
## Warning in att_gt(yname = "taxa_overdose", tname = "year", idname = "Fips", :
## Not returning pre-test Wald statistic due to singular covariance matrix
aggte(overdose_death, type = "dynamic")
## 
## Call:
## aggte(MP = overdose_death, type = "dynamic")
## 
## Reference: Callaway, Brantly and Pedro H.C. Sant'Anna.  "Difference-in-Differences with Multiple Time Periods." Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. <https://doi.org/10.1016/j.jeconom.2020.12.001>, <https://arxiv.org/abs/1803.09015> 
## 
## 
## Overall summary of ATT's based on event-study/dynamic aggregation:  
##      ATT    Std. Error     [ 95%  Conf. Int.]  
##  -0.0288        0.0055    -0.0395     -0.0181 *
## 
## 
## Dynamic Effects:
##  Event time Estimate Std. Error [95% Simult.  Conf. Band]  
##         -10  -0.0015     0.0025       -0.0077      0.0048  
##          -9  -0.0092     0.0067       -0.0260      0.0075  
##          -8  -0.0004     0.0086       -0.0221      0.0212  
##          -7  -0.0002     0.0052       -0.0131      0.0128  
##          -6   0.0080     0.0253       -0.0554      0.0714  
##          -5   0.0104     0.0038        0.0008      0.0201 *
##          -4  -0.0227     0.0222       -0.0784      0.0331  
##          -3   0.0050     0.0066       -0.0115      0.0215  
##          -2  -0.0093     0.0145       -0.0457      0.0272  
##          -1  -0.0096     0.0122       -0.0403      0.0211  
##           0   0.0025     0.0049       -0.0099      0.0149  
##           1  -0.0126     0.0037       -0.0219     -0.0033 *
##           2  -0.0327     0.0098       -0.0573     -0.0081 *
##           3  -0.0351     0.0074       -0.0538     -0.0164 *
##           4  -0.0660     0.0101       -0.0914     -0.0405 *
## ---
## Signif. codes: `*' confidence band does not cover 0
## 
## Control Group:  Never Treated,  Anticipation Periods:  0
## Estimation Method:  Doubly Robust
#gráficos
ggdid(aggte(overdose_death, type = "dynamic"))

3.2 Regressão de uso de maconha

weed <- att_gt(yname = "taxa_weed",
                   tname = "year",
                   idname = "Fips",
                   gname = "Tratado",
                   data = base_2)
## Warning in pre_process_did(yname = yname, tname = tname, idname = idname, : Be aware that there are some small groups in your dataset.
##   Check groups: 2012,2014.
## Warning in att_gt(yname = "taxa_weed", tname = "year", idname = "Fips", : Not
## returning pre-test Wald statistic due to singular covariance matrix
aggte(weed, type = "dynamic")
## 
## Call:
## aggte(MP = weed, type = "dynamic")
## 
## Reference: Callaway, Brantly and Pedro H.C. Sant'Anna.  "Difference-in-Differences with Multiple Time Periods." Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. <https://doi.org/10.1016/j.jeconom.2020.12.001>, <https://arxiv.org/abs/1803.09015> 
## 
## 
## Overall summary of ATT's based on event-study/dynamic aggregation:  
##    ATT    Std. Error     [ 95%  Conf. Int.] 
##  6.493        6.7584    -6.7533     19.7393 
## 
## 
## Dynamic Effects:
##  Event time Estimate Std. Error [95% Simult.  Conf. Band]  
##         -10  -4.7907     2.8072      -12.9850      3.4037  
##          -9  10.9651     3.6991        0.1675     21.7628 *
##          -8 -11.8953     4.8048      -25.9207      2.1300  
##          -7  -5.3488     5.7925      -22.2571     11.5595  
##          -6   9.2674     1.4211        5.1194     13.4155 *
##          -5   4.8605     3.5820       -5.5955     15.3164  
##          -4   8.0872     1.3330        4.1962     11.9783 *
##          -3   0.7674     3.3248       -8.9377     10.4726  
##          -2  13.7616     3.4050        3.8225     23.7007 *
##          -1   1.9593     8.1855      -21.9343     25.8529  
##           0  -8.6802     3.9383      -20.1763      2.8158  
##           1  -1.8370    14.1333      -43.0920     39.4181  
##           2   9.0236     9.7359      -19.3955     37.4427  
##           3  14.6732    14.2983      -27.0636     56.4099  
##           4  19.2854     9.9434       -9.7395     48.3103  
## ---
## Signif. codes: `*' confidence band does not cover 0
## 
## Control Group:  Never Treated,  Anticipation Periods:  0
## Estimation Method:  Doubly Robust
#gráficos
ggdid(aggte(weed, type = "dynamic"))

3.3 Regressão de uso de cocaína

cocaine <- att_gt(yname = "taxa_cocain",
                   tname = "year",
                   idname = "Fips",
                   gname = "Tratado",
                   data = base_2)
## Warning in pre_process_did(yname = yname, tname = tname, idname = idname, : Be aware that there are some small groups in your dataset.
##   Check groups: 2012,2014.
## Warning in att_gt(yname = "taxa_cocain", tname = "year", idname = "Fips", : Not
## returning pre-test Wald statistic due to singular covariance matrix
aggte(cocaine, type = "dynamic")
## 
## Call:
## aggte(MP = cocaine, type = "dynamic")
## 
## Reference: Callaway, Brantly and Pedro H.C. Sant'Anna.  "Difference-in-Differences with Multiple Time Periods." Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. <https://doi.org/10.1016/j.jeconom.2020.12.001>, <https://arxiv.org/abs/1803.09015> 
## 
## 
## Overall summary of ATT's based on event-study/dynamic aggregation:  
##     ATT    Std. Error     [ 95%  Conf. Int.] 
##  0.0943        0.9362    -1.7406      1.9292 
## 
## 
## Dynamic Effects:
##  Event time Estimate Std. Error [95% Simult.  Conf. Band]  
##         -10   0.1968     1.1379       -2.4435      2.8371  
##          -9   0.2638     0.8841       -1.7875      2.3151  
##          -8  -1.3440     1.4136       -4.6239      1.9359  
##          -7   0.9081     2.6160       -5.1617      6.9780  
##          -6   3.7719     2.5132       -2.0594      9.6033  
##          -5   0.4120     1.5634       -3.2155      4.0394  
##          -4  -0.4363     2.8796       -7.1176      6.2450  
##          -3  -0.2460     1.4337       -3.5726      3.0806  
##          -2  -0.3545     1.1041       -2.9164      2.2074  
##          -1  -1.6069     0.7456       -3.3368      0.1230  
##           0  -0.5664     1.2559       -3.4805      2.3477  
##           1   1.2170     1.5075       -2.2808      4.7149  
##           2   2.6266     1.3369       -0.4754      5.7286  
##           3  -0.3460     1.4724       -3.7623      3.0704  
##           4  -2.4598     0.4880       -3.5920     -1.3276 *
## ---
## Signif. codes: `*' confidence band does not cover 0
## 
## Control Group:  Never Treated,  Anticipation Periods:  0
## Estimation Method:  Doubly Robust
#gráficos
ggdid(aggte(cocaine, type = "dynamic"))